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Abstract. We present a method for solving the two- 
dimensional linearized collisionless Boltzmann equation using 
Fourier expansion along the orbits. It resembles very much so- 
lutions present in the literature, but it differs by the fact that 
everything is performed in coordinate space instead of using 
action-angle variables. We show that this approach, though 
less elegant, is both feasible and straightforward. 

This approach is then incorporated in a matrix method 
in order to calculate self-consistent modes, using a set of 
potential-density pairs which is obtained numerically. We in- 
vestigated the stability of some unperturbed disks having an 
almost flat rotation curve, an exponential disk and a non-zero 
velocity dispersion. The influence of the velocity dispersion, 
halo mass and anisotropy on the stability is further discussed. 



scribed in the literature (e.g. Kalnajs, 1977, Hunter, 1992), but 
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The study of the perturbations in stellar disks (spirals and 
bars) has advanced along essentially 2 quite different avenues 
since the early 60 's. 

Numerical N -body simu lations probably offer the most 



flexible tools (e.g. Hohl, 1971; Athanassoula & Sellwood, 1986), 
they are relatively easily set up regardless the complexity of 
the initial conditions and include the description of nonlinear 
evolution. On the other hand, these simulations are very time- 
consuming and suffer from important, but hard to quantify, 
numerical noise. 

A linearized self-consistent mode analysis, which explains 
linear instability during the onset, was the first method used, 
because initially N-body simulations were still essentially in- 
feasible. This approach is somewhat complementary to N- 
body simulations. It heavily relies on simplifying assumptions 
but produces high quality results, of course only in the re- 
gions where the assumptions are valid. Alt hough the litera- 
ture offers general methods (Kalnajs, 1977) to calculate lin- 



ear instabilities, these methods apparently are difficult to ap- 
ply in practical situations. For this reason, several researchers 
have adopted further simplifications in order to handle the 
equations more conveniently, s uch as using cold disks with 
s oftened gravity ( ro^ornre^^ljSl ) and a gaseous approximation 
(Bertin, Lin, Lowe & Thurstans, 1989). 



In this paper we develop an alternative method (section 
2, 3 and 4), which resembles very much other approaches de- 



differs on a few important points. As in many other cases, the 
linearized Boltzmann equation is solved by fourier expansion 
of the perturbation along the unperturbed orbits. However, we 
perform all calculations in coordinate space instead of writing 
everything in action-angle variables. Despite the fact that the 
equations are less elegant, this strategy offers the advantage 
that the perturbed mass density as well as the full perturbed 
distribution is obtained in ordinary space and velocity coordi- 
nates. This considerably simplifies the solution of the Poisson 
equation, since any complete set of potential-density pairs is 
now sufficient to expand the perturbing potential. This is in 
contrast to the action-angle implementations, where usually 
bi-orthogonal sets enter the calculations (see Kalnajs (1977) 
for more details). In addition, the equations can be cast in 
such a form that all calculations are very efficient, since the 
evaluation of the response mass density is essentially reduced 
to the calculation of a Hilbert transform. 

We applied this method to discuss the instabilities of a 
family of unperturbed disk models (described in section 1). 
All models feature a reasonably flat rotation curve and an ex- 
ponential disk. This family consists of almost isotropic distri- 
bution functions, with varying velocity dispersions. Stability 
increases with increasing velocity dispersions and increasing 
halo mass (section 5). Section 6 sums up. 

1. The models 

Before presenting the method of the mode calculation in more 
detail, we first introduce the unperturbed galaxy models for 
which the calculations were performed. All our models have 
the same unperturbed potential and mass density, but feature 
different distribution functions, with varying velocity disper- 
sion, streaming velocity and anisotropy. 

1.1. The unperturbed potential 

The potential of the unperturbed galaxy was constructed as a 
sum of two Kuzmin-Toomre disk potentials with different core 
radii. In the plane of the disk, it is given by 
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Note that potentials are defined as binding energies, with 
a positive sign. This potential produces a rotation curve which 
is much flatter than a single component Kuzmin-Toomre po- 
tential (see fig. |l|). The ratio of the flat part to the rising part 
is about 6/1. 
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Fig. 1. Rotation curve of the unperturbed galaxy potential 

Although the rotation curve tends to rise somewhat too 
slowly near the centre, it has, at least qualitatively, a realistic 
behaviour. 
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Fig. 2. Total mass inside r of the disk (full curves) and the halo 
(dashed curves) for H/ D=2.5, 5 and +oo. 



1.2. The mass density profile 

It is well-known that strongly flattened galaxies usually have a 
substantial amount of dark halo mass, extending much further 
than the visible component. It is sufficient for our purpose to 
model it by a spherical and pressure-supported galaxy compo- 
nent. Therefore it is reasonable to assume that this halo com- 
ponent only influences the stability behaviour of the disk by its 
contribution to the global potential. The same roughly holds 
for the central bulge, which is hot and has a three-dimensional 
structure as well. Thus both the halo and the bulge are "inert" 
and the corresponding potential and mass density are taken to 
be spherical and denoted by Vo,h(t) and po,ii{r). 

The disk itself is the only component which is supposed to 
consist of "active mass", sensitive to instabilities. The unper- 
turbed disk is supposed to be two-dimensional and axisymmet- 
ric, with a surface density po,D and a potential Vo,d- We chose 
an exponential mass profile with a central core: 



Of course, the halo mass density should everywhere be pos- 
itive. This puts an upper limit on the value of a in (|^). The 
relative contribution of the disk and halo components are quan- 
tified using the halo-to-disk factor, H/D, which gives the pro- 
portion of the total mass inside the radius r max (=6) for the 
halo and the active disk. Self-consistency with a non-negative 
spherical halo puts a lower limit on H/D of about 2.5. Note 
that the assumption of a spherical halo is not a crucial one. If 
one would assume an oblate halo, the only effect would be a 
somewhat lower minimum H/D. 

1.3. The distribution function 

We will examine the stability behaviour for different stellar 
distributions. According to Jeans' theorem, the unperturbed 
part of the distribution is a function of two integrals of motion, 
the binding energy E and the angular momentum J, defined 
by 
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Although this mass density extends up to infinity, the ac- 
tual models have an outer limit at r = 6, and the mass density 
reaches zero at that point. This we achieve by fitting the distri- 
butions (^) with several finite components. The relative error 
made at the outer edge is only of the order of 0.1%. The pa- 
rameter a determines the total disk mass. 

Since the total system should be self-consistent, the total 
potential in the plane of the galaxy is given by the sum of the 
disk and halo component: 



Vo(r) = V , D (r) + V , H (r). 



(3) 



Since the total potential Vo has a fixed form (|l|) and the 
potential of the disk is determined by the surface mass density 
(^), we can calculate Vq,h = Vo — Vo,d- Thus follows the halo 
mass density 
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In order to generate a variety of finite disks, based on the 
potential (|l]), the distribution function is written as a linear 
combination of basic distributions: 
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All components only have a (everywhere positive) contri- 
bution for orbits lying completely inside r max , i.e. for the region 
where (see also fig hq) 
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with r max the radius of the edge of the disk. Equation (g) is 
required in order to exclude the region r_ > r m ax. In addition, 
the distribution goes to zero at this edge in a smooth way, so 
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that the first derivative remains finite everywhere. The expiicit 
form of the components is fisted in the Appendix. 

The expansion coefficients ct are determined by a feast 
square fit of the corresponding mass density to the proposed 
exponentiai form (Q) (the coefficients are forced to be positive, 
in order to avoid negative distributions). By choosing an ap- 
propriate set of components /o,t, we were able to create the 
desired orbital densities. In all the models, the error on the fit 
to the mass density never exceeds 1% of the centraf vafue. 

We constructed 4 modefs, fabefed I to IV (the expficit form 
of the distribution functions is fisted in the appendix). Along 
the sequence, the models become more and more rotation- 
supported, having an increasing streaming velocity and de- 
creasing temperature. Fig. ^| shows the streaming velocities 
and dispersions for the cofdest and hottest case. The disper- 
sions ail go to zero at the edge of the disk. Note that model I is 
perfectly isotropic and has a iineariy increasing mean veiocity 
curve. For all disks, To omre's local axisymmetric instability 
criterion (Toomre, 1964) 



Q = 
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(with o r the radial velocity dispersion and k the epicyclic 
frequency) is everywhere higher than 1. 
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which leads to 
= r 2 + V (r+) ~ r 2 _V {r-) 

and 
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Inversely, r + and r_ are found as the roots for r of (JlfJ) 
for a particular E and J. We preferred these variabies over the 
normal [E, J) space, not only because they have an easy physi- 
cal interpretation, but also because this representation is more 
related to our method for solving the linearized Boltzmann 
equation, which employs a grid interpolation in the turning 
point space (see foffowing section). 

We have chosen finite disks in order to compactify phase 
space. However, as can be seen from the distribution function 
(fig. ^), the disk reaches this hmit at r max in a very smooth 
way. This is important since it has been proven that a sharp 
edge or, more generaily, a sharp feature in (E, J) space, can 
introduce additional in s tabilities which might n ot always be 



physical (Toomre, 1964, sellwood & Kahn, 1991) 




Radius 

Fig. 3. Streaming velocity (full line), azimuthal velocity dispersion 
(long dashed line) and radial velocity dispersion (short dashed line) 
for Model IV (top panel) and Model I (bottom panel). 

In fig. ^, the distribution function of modef III is shown 
in turning point space. The turning points of an orbit are the 
largest (resp. smallest) distance from the centre that an or- 
bit can reach, and are caifed apocentre r+, or pericentre r_. 
By convention, we take r + always positive, while r_ has the 
same sign as J. Evidently, these quantities are integrals of the 
motion and r+ > r_j (on circular orbits, r+ — r_). For a 
given pair (r+,r_), the energy and angular momentum follow 
immediateiy from 



E = V (r+, 
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Fig. 4. Turning point representation of the distribution function of 
model III. 



2. The Poisson equation 

Since we want to study uniformly rotating perturbation modes, 
we suppose the following form for the perturbing potentiai: 



V'(r,9;t) = V'(r)e 



i(m6 — u)t) 



(14) 



In this expression (and in all later similar expressions), only 
the real part corresponds to the physical quantity. There is no 
problem though to calculate with these complex expressions, 
since all equations are linear. This perturbation is rotating with 
a pattern speed of Re(w)/m and is exponentially growing with 
a growth factor \m(ui). A more general pattern can be obtained 
as a superposition of such components with different m. 

We wiil search for the linear modes using a matrix method, 
which requires a set of potentiai-density pairs suitabie for the 
expansion of the perturbation. The liter ature offers a variety 
of anatytica t sets for infinite disks (eg. |Ctutton-Brock, 1972 , 



(11) 



Qian, 1993) as welf as for finite ones (eg. Hunter, 1963). The 



choice of a particuiar set largely depends on the structure of 
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the unperturbed potential Vo and density po- The importance 
of a suitable set of potential-densities is clearly illustrated by 
the displacement mode of a self-consistent disk. If the model 
has no inert halo component, a simple translation of the system 
is a valid "perturbation" . In linear theory, th is cor responds to 
&n m — 1 mode of the form (see also section 4.1.S) 



(15) 



The potential-de nsity pairs shou ld of course be able to fit 
this behaviour well (Weinberg, 1991). Moreover, it is very im- 
portant to choose the set as efficiently as possible, so that not 
too many expansion terms need to be calculated. For these 
reasons, we decided not to rely on existing analytical sets, but 
to create numerical sets which are tailor-made for the present 
unperturbed disk. If the maximum radius of the finite disk is 
denoted by r max , the set of densities for m = is given by 



r 

^max 



Ps.nW = Po,r>(r) * cos(7rn- ), 

and for ra > by 



m / \ r, 

psA r ) = r 



-i dpo.D , , , r 
-(r) * cos(7rn 



dr 



(16) 



(17) 



These densities have the right behaviour in the centre and 
at the edge of the galaxy. In addition, their number of radial 
nodes is equal to n. The corresponding potentials V<^ n are eas- 
ily calculated numerically using the integral form of the Poisson 
equation. These potential-density pairs prove to be very effi- 
cient for the expansion of the normal modes. Of course, these 
modes in principle could have been expanded using analytica l 
complete sets described in the literature (e.g. Hunter, 1963 ), 
but this would require a long summation up to very high or- 
der, particularly because the density falls to zero very smoothly 
at the edge. 



3. The linearized Boltzmann equation 

3.1. Integral form 

The total disk is modeled using two parts, a time-independent 
axisymmetric unperturbed part, and a small perturbation. 
Hence, the distribution function is written as 



f(r, 8, v r , v g ; t) = f (E, J) + f'(r, 6, v r ,v e ;t), 



(18) 



while the total potential, which is defined as a binding en- 
ergy (being positive everywhere), reads 



V(r,e;t) = V (r) + V'(r,9;t). 



(19) 



If the system is exposed to the perturbing potential V', 
t he corresponding linearized response distribution is given by 
(Vautcrin fc Dejonghe, 1995) 



ri _ dfo ( , dfo t , 

where f' E and f'j satisfy (using Poisson brackets) 
^-[f' E ,E] = [E,V'\ 



and 



dt 



J --[f'j,E] = [J,V']. 



(20) 



(21) 



(22) 



Integration with respect to the time along unperturbed or- 
bits of both equations and expansion of the Poisson brackets 
immediately yields (since the left hand side is the total time 
derivative along an unperturbed orbit) 



f B {r ,0 ,v T ,v e ;0) = \ {v r — + — -^-)dt 



/S(r°,0> r °,^;O) 



dV 



dt. 



(23) 
(24) 
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The integrand is to be evaluated along the unperturbed orbit 
which passes at t = through the point (r°, 8°, v®, Vg). One 
should note that these expressions can only hold if the per- 
turbation was zero at t — — oo and is, at least infinitesimally, 
growing. 

With a potential of the form 



V'(r,6,t) = V'(r)e 



we have 



(25) 



f I dV' V' . i( m 8-ut) ,„„>. 

IE = I ( V r-^T + imv 0—) e dt (26) 



f'j = / imV'e^-^dt. 



(27) 



These integrals are to be evaluated along unperturbed orbits, 
so one needs a way to handle such orbits. 

3.2. Fourier expansion along unperturbed orbits 

Since the unperturbed potential Vo is axisymmetric, the cor- 
responding orbital structure is particularly simple. Taking ad- 
vantage of the conservation of angular momentum J, one can 
conclude that the radial coordinate behaves like the coordinate 
of a particle in a one-dimensional (so-called effective) potential 



(28) 



Supposed that the orbit is bound (which is the only case 
of interest for galaxies), the radial coordinate along an unper- 
turbed orbit should be a periodic function of time, with angular 
frequency uj t . From this it follows immediately that v r and vg 
are periodic with frequency uj r as well. As the mean value of 
vg should not necessarily be zero, the angular coordinate 9 will 
be a superposition of a periodic function and a uniform "drift" 
velocity: 



L0gt + 9 p {t) 



(29) 



(with 9 P a periodic function with period uj r ) . 
We will now rewrite the integrands of (^6|) and (^i|) by 
factorizing the part which is periodical with angular frequency 



f' E = 
f'j = 
with 
I E (t) - 



iBity^o-^dt 



Ij(t)e 



i(mu 9 -u)( 



dt 



v r{t)-^r{t) + tmv e {t)^-\e 



and 

Ij{t) = imV'{t)e mB * {t) . 



(30) 
(31) 

(32) 
(33) 



Since I E (t) and Ij{t) are periodic, they can be expanded 
in Fourier series: 



I E (t) = ]T I E ,ie 
i=— W 

imax 



j,ie 



(34) 
(35) 



When these forms are substituted in ( p0| ) and ( p0| ) , the integra- 
tions can be carried out analytically, at least if the perturba- 
tion is growing, Im(a>) > 0. Note that, when the perturbation 
is decaying in time (Im(oj) < 0), one can still obtain a solution 
by performing the integrations ( |30| ) and ( ^o| ) from t = to 
t = +oo. 



3.3. Grid interpolation 

The Fourier expansion strategy yields the response distribu- 
tion for any given perturbing potential, but a number of steps 
are involved even for calculating one single point of the distri- 
bution: the orbit should be integrated for at least one half of 
the radial period, and the appropriate functions are to be ex- 
panded in Fourier series. Since the perturbed distribution will 
be evaluated many times (e.g. for the calculation of the per- 
turbed mass density one has to integrate over the velocities), 
it is absolutely necessary to be able to evaluate the response 
as fast as possible. Our solution is to calculate all appropriate 
parameters (the frequencies ui r and uig and the Fourier expan- 
sions Ie,i and Ij t i) on a (two-dimensional) grid in integral space 
and to store the result on disk for future use. We found that 
is was not so convenient to use E, J as grid coordinates, since 
this grid is very inhomogeneous and has a relatively low den- 
sity near the circular orbits, where disk galaxies have a dense 
population. Moreover, the circular orbit limit on this grid is 
usually not given in an analytic form. For these reasons, we 
used r+ and x = r_/r+ as grid coordinates. In this represen- 
tation, we can apply a simple rectangular grid which is very 
dense close to circular orbits and in the neighbourhood of the 
centre of the galaxy. 

The (r+, x) space is, up to the maximum radius r+ = r max , 
discretized using a rectangular grid. On each point of this grid, 
the Fourier expansion for the orbit starting at t = in its 
apocentre (by choice) is calculated and stored in a library, to- 
gether with a tabulation of t[r] and 9 p [r]. In the following para- 
graphs, we will frequently use both the time and r-dependence 
of the same variable. In order to avoid confusing notations, in- 
dependences will be written using square brackets. The time 
and coordinate system is chosen in such a way that 9 and t are 
zero at the apocentre. 

When the response distribution is later required in a point 
p° — (r°, 9° — 0, , Vg; t° = 0) in phase space, the correspond- 
ing turning points r° and r°_ are calculated and all parameters 
are interpolated using the four closest points on the grid (note 
that, since the response is a known periodic function of 9, it 
is sufficient to evaluate it at 9 = 0). In this way, we get an or- 
bit out of the library with the correct integrals of motion, but 
passing through its apocentre at t = 0. We will denote this li- 
brary orbit and all corresponding quantities with a superindex 
L. 
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Fig. 6. Interpolation grids of the lowest order Fourier coefficients 
for a Vg j perturbation (only the x > part is shown) 



As the point p° is in general not the apocentre of the orbit, 
one should take in account an offset in time for the actual orbit 
passing through p°: 



lk(t + t L \r 



(36) 



Note that t [r ] is actually a double-valued relation. How- 
ever, the sign of u° can be used to determine which branch of 
the relation should be used. In general, there is also an offset 
in azimuth: 



0(t)=wft + #(t + t>°])-#[r o ], 



(37) 



which is of course the reason why the functions t[r] and 
6 P [r] are stored. Note that the presence of these functions does 
not reduce the speed of the calculations, but implies only a 
slight increase in used disk space. The resulting perturbed dis- 
tribution in that point follows then immediately: 



E 



J(l^t L [r o ]-m0^[r }) 



i(liUr + mujg — lo) 



(38) 



and 



f'j = 



(max j-l e i{l^t L [r°]- m e^[r ]) 
> — . 



(39) 



The Fourier coefficients which are stored on the grid only 
depend on the unperturbed and the perturbing potential. Fur- 
thermore, since the latter will be expanded in a basis, this has 
to be done only for the discrete set Vg^. The orbits are ob- 
tained by integrating the system over half a period using a 
fourth order Runge-Kutta. Special attention should be paid to 



r- = (or J = 0) orbits, which have to be integrated in rectan- 
gular coordinates. In addition, sinc e the Fourier c oefficients are 
discontinuous at r_ = (see also Kalnajs, 1977), this should 



be performed twice, for "infinitesimally small" opposite values 
of J in order to get the correct left and right limit. 

The Fourier expansion is performed up to the order Z max = 
14. This number should be relatively large in order to provide 
an accurate fit to the high excentricity orbits, which show fast 
variations close to the centre. A grid resolution of 60 x 60 proved 
to provide more than enough accuracy, particularly because the 
functions do not show very steep gradients in the turning point 
space. In these circumstances and for a set of 8 perturbing 
potentials, the calculations take less than 4 hours on a 100 
Specfp92 workstation. This work only has to be redone when 
one changes the unperturbed potential. 

Fig. M shows some of the Fourier coefficient grids for a m — 
2 perturbing potential. It clearly illustrates that the higher 
order Fourier terms are only important for small x — r-/r+ 
and for large r+, which means for very eccentric orbits. This 
indicates that, for relatively cool galaxies, a small number of 
Fourier terms can already yield reliable results. 

Fig. |?] displays some response distribution functions for dif- 
ferent uj. The left and middle panels have the same pattern 
speed Re(o;)/m and are taken somewhat inside the corotation 
radius. The left panel shows a response with a much smaller 
growth rate and therefore has a more prominent corotation res- 
onance, appearing as a ring-shaped feature. This picture still 
has to be multiplied with dfo/dE which, for relatively cool 
disks, will emphasize the stars in the neighbourhood of circu- 
lar orbits. One of the most crucial effects of non-zero velocity 
dispersion is that the distribution function populates stars on 
all possible orbits that are resonant with the perturbations, 
not only the resonances on circular orbits. These orbits fill a 
complicated region in phase space, cuts of which can be seen 
in velocity space in fig. ^. Due to the differential rotation, a 
substantial number of orbits at this radius are in corotation 
resonance, although this is not the corotation radius. 




Fig. 8. The resonant values of u at r = 1 in the (vg,v r ) plane for 
the three most important m = 2 resonances. Circular velocities are 
indicated. 

Fig. H shows, for a fixed radius, the velocity dependence of 




Fig. 7. The real part of the perturbed distribution Re(f' E ) at r° = 1.8 for a typical m = 2 perturbation. From left to right: at ui = (0.7, 0.07), 
u! = (0.7,0.7) and u> = (1.5,0.1) (dotted regions are negative responses; the small ticks on the vg axis correspond to circular velocities). 



the resonant w value for the dominant m = 2 resonances. This 
value is given by 



to = lw r + mujg. 



(40) 



The lowest order resonances for stars rotating in the direct 
sense (wg > 0) are 1= -1 (ILR), (CR) and 1 (OLR). How- 
ever, there is a problem in this definition for counterrotating 
stars, since wg discontinuously changes its sign at vg — 0. This 
can be resolved by noticingthat if, for counterrotating stars, 
I is replaced by I + m in (ftOh, the resulting w has an overall 
smooth behaviour. This is a consequence of the fact that for 
zero momentum stars w r — 2wg. The ring-shaped CR position, 
as shown in fig. ^ (for varying w), can easily be found back in 
the left panel of fig. (for a specific value of w) . 

Since w r and wg depend mostly on the energy E, the 
CR, OLR and higher resonances are roughly concentric circles 
around the origin. At the ILR though, the sum largely cancels 
out and result is mostly dominated by the angular momentum 
J. 

3.4- Integration over the velocities 

The method of the previous section can be used to calculate 
the perturbed distribution in a particular point of phase space 
(r, 9 = 0, v r ,vg). Writing the velocities (v r , vg) in polar coordi- 
nates (v,a), this results in a response 



f'(r,e = 0,v,a,t = 0) = 



with 

Pi (r, v, a) = lw r + rawg . 



E 

i=-ima 



Ai(r, v, a) 
pi(r,v,a) - w' 



(41) 



(42) 



In order to obtain the perturbed mass density at a radius r 
and for 8 — 0, one has to integrate this expression over the ve- 
locities up to the escape velocity. Since the position of the pole 
Pi in general depends also on the velocities, the denominator 
should be kept inside this integration: 



p'(r) = 



i=-!ma 



Mr 



Pi(r,v, a) - w 



vdvda. 



(43) 



In each term of the sum, we will switch the integration vari- 
ables from (v,a) to (pi(r, v, a), a). The velocity then becomes 
a function v = Vi(r,pi,a) and the integral is written as 



o'(r)= £ 



ima 

l 

i = -imax"' P! . 



Pi, 



Pi — u 



Ol X 



dv 



Ai(r,Vi(r,pi,a),a)Vi(r,pi,a) — (r,pi,a)da. (44) 
o °P' 



The limits of the first integral can be chosen freely, as long as 
the interval is large enough to cover all poles. If pi (r,v, a) is not 
a monotone function of v over the entire region of interest, this 
transformation can still be performed by splitting up the inte- 
gral in parts. There is a simple reason to write the expression 
in this way: for a constant r, the inner integrals are functions 
which only depend on pi. If we sum all these functions over I 
and denote this summation with W(r,p), the perturbed mass 
density reduces to 



p(r) = 



W(r,p) 



(45) 



which is a simple weighted integral over all real pole posi- 
tions. This integral represents the Hilbert transform of W with 
respect to p, assuming that p is also a function of w. For each 
value of r, the function W(r,p) can be tabulated and stored 
for later use. The evaluation of the perturbed density for any 
value of w is then reduced to the numerical evaluation of a 
one-dimensional integral, which can be performed very fast. 

In practice, W(r,p) is not calculated using the explicit inte- 
gral d44j), because this would require the (piecewise) inversion 
of pi(v), which can be awkward. Instead, an alternative method 
is used, based on the two-dimensional integral (^3|). For any 
value of r, W(p) is maintained in a tabular form, containing 
function values at fixed points with interleave Ap. 

A numerical integration of the two-dimensional integral 
( ff3| ) is performed, using a rectangular grid over [0, v csc ] x [0, 27r], 
with interleaves Ad and Aa. In each cell of this grid, the in- 
tegrand of ( |43[ ) is approximately constant with respect to v 
and a. Therefore, the contribution of this cell to the integral 
is a pole function with respect to w, having a pole position 
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nances, particularly for the corotation (CR) and outer Lind- 
blad resonance (OLR). The inner Lindblad resonance (ILR) 
is, even for hot disks, much sharper, and remains roughly at 
t he same value fo r u) throughout the galaxy. As noticed earlier 
(Lin et al., 1969), this interesting behaviour follows from the 
structure of the unperturbed potential. Due to this behaviour, 
the ILR plays a crucial role in the stability behaviour, even for 
high velocity dispersions. 

Note that, for the panel corresponding with r — 1, one can 
easily correlate the regions in p spanned by the contributions 
of the individual resonances with the variation in ui shown by 
the corresponding contour plot in fig. pi 



4. Construction of normal modes 

A matrix meth od is applied for searching the normal modes 
(Kalnajs, 1977). A general perturbing potential is written as 



V 



E 

i=0 



V-r> 
- s 



(46) 



Of course s should be large enough so that all desired modes 
can be expanded accurately (the limitation on s puts a limita- 
tion on the oscillatory behaviour of the modes). 

Using the methods described in the previous paragraphs, 
the pole density functions Wg\(r,p) associated with the per- 
turbing potentials Vcf\ are calculated and expanded as 



E 

3=0 



C 3,i(P)P 



S,3- 



(47) 



Fig. 9. The pole distribution W(lu) for the perturbed mass density 
using model III and a V s 2 perturbation. The bumps at negative p 
correspond to a counterrotating resonance. 



P = Pi(r,v,a) and an amplitude / = Ai(r, v, a)vAvAa. This 
amplitude is now added to the table value of W(p) that has 
a p value lying closest to P. When all cells are summed, di- 
viding all values in the W(p) table by Ap yields a numerical 
tabulation of W(p). In practice, the summation over the cells 
is refined using a trapezoidal rule. 

The values of the function W(r, p) are stored on a grid 
containing 40 points in the r dimension and 10 4 points in the 
p dimension, and the integral (^) is calculated by linear in- 
terpolation between the points in the p dimension. However, 
one should be careful with values of uj lying in the neighbour- 
hood of the real axis. If Im(u) is zero, the result is meaningless 
when W(r,p) 7^ in a region [Re(p) — e, Re(p) +e] (in this case, 
stars are in resonance with a steady perturbation and the lin- 
ear approximation breaks down anyhow). In addition, for such 
a region of non-zero weight for poles, |Im(w)| should be large 
in comparison with the grid distance of the p grid in order to 
have a reliable interpolation. This is the reason why we take 
such a dense grid for p. 

In fig. [| the pole distribution W(r,p) at various radii for 
a typical configuration is shown. For small r, the contribu- 
tions of each resonance are easily distinguished, while for large 
r the behaviour becomes more complex since the stars move 
slowly and many higher order resonances occur. This picture 
also clearly illustrates the non-localized character of the reso- 



This expansion is obtained using a least square fit. One can 
easily obtain c™j(p) in tabular form by performing the fit on 
each of the rows in the r direction of the W(r, p) grid. We can 
now write the response density p' s t for each of the potentials 

V s ™as 

3 

p' 7 s,i = 2Z c TA u )p%h ( 48 ) 
3=0 

if we define 

# t ( W ) = / -^dp. (49) 
The response potential of ( |46[ ) now follows immediately: 

E (E r ) r: :- ( 5 °) 

j=0 \i=0 J 

For a self-consistent mode, this response should be equal to 
the original perturbing potential (|46[). So the search for normal 
modes is reduced to the search for those values of u) for which 
the matrix C(uS) (with elements c g , p ) has a unity eigenvalue 
A. The corresponding eigenvector gives the expansion of the 
normal mode. 

A typical situation for the multi- valued eigenvalue relation 
\{lo) is shown in fig. The steep edge in Re(A) and the 
sharp peak in Im(A) at small positive pattern speeds are the 
dominating features and are caused by the ILR because the 
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Fig. 10. Real and imaginary part of the multi-valued eigenvalue 
relation at Im(u;) = 0.11 for the model III. The different kinds 

of lines correspond to different branches of the relation. The unique 
eigenvalue with A = 1 is indicated. 



reaction of the disk is strongest at ILR (see also fig. |^). For 
each value of Im(u;) and for each branch, there appears to be 
only one (small and positive) value of Re(o;) where Im(A) is 
zero. The value of Re(A), which decreases for increasing Im(u;), 
further determines the mode location. 

Obviously this search for unity eigenvalues in the complex 
plane has to be performed numerically. This is the main reason 
why we put so much emphasis on the development of a method 
to evaluate C(u>) as fast as possible. It takes less than a second 
to calculate the eigenvalues for a particular uj, using an expan- 
sion up to an order a = 8. Due to the efficient structure of the 
potential-density pairs, this limit is certainly sufficient for an 
accurate expansion of the unstable modes, since they prove to 
be of a relatively low order. The program calculates a detailed 
numerical tabulation of the complete dispersion relation in a 
sufficiently large region of the complex plane. This tabulation 
is further used to determine all modes using a binary inter- 
section algorithm. This whole procedure takes less than one 
hour. 

4-1. Checkpoints 

Of course, there is a strong need for good checkpoints before 
starting to interpret the results coming out of this method. 
We checked the system in three ways, together covering more 
or less the whole chain of operations. The positive results of 
these checks are not only a good indication that the calcula- 
tions are not obviously wrong, but they provide a good tool 
to determine the values of the various parameters (grid resolu- 
tions, expansion limits, ...) in order to get results with sufficient 
accuracy. Most of these parameters were actually determined 
experimentally, using these test cases. 



4.1.1. Comparison with direct integration 

For a given perturbing potential, the linearized Boltzmann 
equation can be integrated numerically, at least for sufficiently 
fast growing perturbations. This yields values for the perturbed 
distribution function which should be the same as those coming 
from the Fourier expansion along the orbits. 

4.1.2. Uniformly rotating disks 



In a previous paper ( Vauterin fc Dejonghe, 1995| ), the mode 



analysis for uniformly rotating disks has already been per- 
formed in an independent way. If the same models are treated 
using the present approach, the same results should come out. 

4.1.3. The displacement mode 

As mentioned already, for models without passive component, 
a simple displac ement of the disk is a valid "perturbation" 
( Weinberg, 1991 ). In a linear approach, the responses corre- 
sponding to this mode are obtained by derivation of the un- 
perturbed values with respect to a rectangular coordinate (e.g. 
x). For the unperturbed mass density, this results in (using the 
fact that x = re 19 ) 



dpo.D , dpo,o , s i0 

— a ( r = — a — v e • 

ox or 



(51) 



The distribution function gives rise to the following pertur- 
bation: 



dfodE dfodJ 
8E dx dJ dx' 

so that we immediately have that 



,1 _ UCi_ _ UVQ j, 

dx dr 



f'j 



and 
dJ_ 

dx 



(52) 



(53) 



(54) 



We used a self-consistent model based on a single Kuzmin- 
Toomre potential to prove that the method is able to find such 
modes. Note that, although this displacement mode occurs at 
lu — 0, there is no problem with the resonances since the pole 
distribution W(r, 0) = for m = 1 perturbations, which we 
verified numerically. 

5. Results and discussion. 

We performed the linear mode analysis for the models I to IV. 
These models have a roughly isotropic distribution function, 
with a velocity dispersion which is decreasing from model I 
to model IV. The aim of this set of models is to address the 
dependence of the stability on two important parameters, the 
velocity dispersion and the halo to disk proportion, H/D. Both 
parameters are now widely believed to have a crucial impact 
on the stability. The velocity dispersion of a particular model 
will be represented by its central value, <T ccnt , which is a good 
estimator of the relative overall behaviour (see fig. ^ . Unstable 
modes with m varying from to 4 were sought in all models 
and for varying H/D. One can easily obtain the results for 
different H/D simply by varying the total magnitude of the 
unperturbed distribution function, represented by a in (^). 
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Fig. 11. Instability limits in the <7 C cnt vs. H/D plane (taken as 
Im(dj) = 0.005) for the lowest order modes in models I to IV. The 
lower left corner corresponds to unstable regions. The models are 
represented by small circles, and the model number decreases from 
the left to the right (the lines between the models are cubic spline 
interpolations). 



Fig. 12. Contour lines of the growth rate of the dominating m = 2 
mode as a function of o" CC nt. The small circles represent the model 
points (the lines connecting the models are interpolations). The dot- 
ted line is an extrapolation at Im(oj) = 0, and is the same as the 
height m = 2 limit in fig. 13. 



In fig. [y], the stability limits for these modes are plotted in 
the (T ccn t vs. H/D plane (we chose Im(aj) = 0.005 as a stability 
limit, since Im(u) = is unreachable with our method). These 
curves were obtained by a cubic spline interpolation between 
the four models, which are represented by small circles. For 
most of the modes, only a small fraction of the models can 
actually reach the unstable regime. The other models only be- 
come unstable if they are combined with a spherical halo which 
is not everywhere positive (and hence H/D < 2.5). The curves 
in figure hi] and [12j are always interpolations between all four 
models, but many models appear under the H/D — 2.5 limit. 

The region under each curve is the unstable region. 
At every point of the figure, the actual (in)stability of 
the disk is of course determined by the highest limit. 
Not surprisingly, the disk is effectively stabilized by in- 
creasing the velocity dispersion or by adding mass to 
the inert halo, or a combination of both. This be- 
haviour has already been r eported for various other cases, 
such as Kuzmin potentials ( Sellwood, & Athanassoula, 1986 



Athanassoula fc Sellwood, 1986|) and quadratic potentials 



( |VaTiterhT7riOcjongh^*^99*5 ) . Another feature, which was also 



shown to be present in quadratic potentials, is that the slope 
of the stability limits increases for increasing m (at least for 
m > 2). For hot disks, the m = 2 mode is the only unstable 
one, while at the cool end, m = 4 turned out to be the most 
unstable perturbation. Note that the model I was found to be 
stable for all calculated harmonics. 

The growth rate for the most unstable m = 2 mode is 
shown in fig. [l2|, again as a function of <T ccnt and H/D. This 
picture clearly shows that the growth rate tends to increase 
in a more than linear way as the halo mass decreases. There 
is also a tendency that the contour lines become less depen- 
dent on the velocity dispersion for large values of the growth 
rate. This can be understood from the fact that varying the 
velocity dispersion re-distributes the weights W(p) of the pole 



positions, and does not change the overall intensity (which is 
influenced by the factor H/D). When Im(oj) becomes large, the 
pole distribution is "seen" from a big distance, and the internal 
positions become less important than the overall intensity. 




1.6 1.8 
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Fig. 13. Growth rate of the m = 2 mode as a function of Toomre's 
Q factor in the centre of the galaxy. 

Athanassoula and Sellwood (1986) suggest that a value of 
2.0-2.5 for Toomre's local stability parameter Q might by a 
useful criterion for stability of the disk a gain st global m — 2 
perturbations. For the present models, fig. O plots the growth 
rate of the dominant m = 2 mode against the corresponding 
Qccnt, the value in the centre of the galaxy. As Q is an in- 
creasing function of the radius, Q ce nt is a lower limit of any 
averaging of Q over the disk. For each modeLthe variation in 
Q is obtained by changing H/D. Equation (hoi) shows that Q 
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simply scales with po, which is influenced by H/D. The fig- 
ure shows that a value between 2.0 and 2.4 is indeed again a 
good estimation of the stability limit. Our values are slightly 
larger than those of Sellwood & Athanassoula, but this can be 
explained by the fact that they used a softened gravity. 
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Fig. 14. Pattern speed of the perturbation as a function of the 
growth rate. The dashed horizontal line corresponds to the limit of 
the presence of an ILR. 

The pattern speed of the perturbation, given by Re(u>) /m, 
is plotted against the growth rate in fig. [14] (note that Im(tj) = 
is, again, an extrapolation). This curve is shown for the 
three unstable models II, III and IV, and the variation in 
the growth rate is obtained by varying H/D. The endpoints 
of these curves are a consequence of the self-consistency re- 
quirement. The pattern speed has a clear increasing correlation 
with the gro wth rate, which is again in agreem ent with N-body 
simulations (Athanassoula & Sellwood, 1986). In addition, this 



figure shows that, for the same growth rate, the pattern speed 
decreases with the velocity dispersion. On the same picture, 
the upper limit on the pattern speed for the presence of an 
ILR is shown as a horizontal line. From this, it is clear that 
none of our models have an ILR, again in agreement with the 
results decribed by Athanassoula & Sellwood (1986). In fact, 
orbits which are in ILR cause such a violent response (see e.g. 
fig. ^| and fig. [lo]) that it would be very hard for a galaxy to 
maintain much of them. 

As a matter of illustration, fig. [l5| shows the density profile 
of the dominating m — 2 mode in model III. This perturbation 
develops its usual behaviour, with a bar-shaped central part 
and losely wound spiral arms in the outer regions. The spiral 
structure is always found to be trailing. 

6. Conclusions. 

A big part of this article was devoted to the description of 
a method for finding linear modes in stellar disks. The pro- 
posed strategy heavily relies on existing techniques, such as 
the matrix method and Fourier expansion along the unper- 
turbed orbit, but differs from previous approaches by the fact 
that everything is calculated in ordinary coordinate and veloc- 
ity space and that a numerical set of potential-density pairs 




Fig. 15. Dominating m = 2 mode for the model III at H/D = 2.5. 
The two circles correspond to the CR and OLR. Both the streaming 
velocity and the pattern speed are clockwise. 



is used. With the proposed scheme, the full perturbed distri- 
bution function is obtained with no extra calculation costs. In 
this way we have shown that calculations in coordinate space, 
although much less suited for theoretical considerations, can 
offer a fast and flexible alternative to action-angle variables 
when it comes to numerical computation of normal modes. 

This method was applied to a set of unperturbed disk mod- 
els, having a more or less realistic potential and an exponential 
mass density. This disks are embedded in a spherical, inert halo 
in order to obtain a self-cons istent model. In agreement with 



various other studies (e.g 



Athanassoula & Sellwood, 1986 



Vauterin fc Dejonghe, 1995), the calculations showed that 



these disks can be stabilised by increasing the velocity dis- 
persion and/or the halo mass. For almost isotropic velocity 
dispersions, a Q ce nt value of 2.0 — 2.5 turned out to be a rea- 
sonable stability limit. 

Comparison of the stability of the present models with 
t he behaviour of disks em bedded in quadratic potentials 



( Vauterin fc Dejonghe, 1995 ) shows a striking resemblance. 
Qualitatively, the stability behaviour of those simple uniformly 
rotating systems shares all the features that we have found for 
the more sophisticated models discussed in this paper. And, 
to a certain degree, there is even a quantitative agreement. 
It seems that, for some important aspects of the stability be- 
haviour, the structure of the unperturbed distribution might be 
more inportant than the nature of the unperturbed potential. 
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A. The unperturbed distribution functions. 

The simple r+ = r ma x limit of the unperturbed distribution in 
phase space, defined in E-J space by 



/in = Et fiM x ( 5.69 x 10" 
+7.69 x 10" 
+1.40 x 10" 



/iv = -El, o.oi x ( 



6.69 x 10" 
+7.90 x 10" 
+7.54 x 10" 
+4.29 x 10" 



£ n e 4J 
£ 8 pw( J, 2) 
£ 6 pw( J, 4) 



E 25 e 9J 

£ 25 pw(J,2) 
£ 20 pw(J,5) 

£ 14 pw(J,7) 



(A8) 



(A9) 



The function pw(i, n) which we introduced equals x n for x > 
and zero for x < 0. If n is integer and larger than one, this 
function has a continuous first derivative. 



E L (E,J) = E + 

ZVrr 



Vo(r max ) > 0. 



(Al) 



is not so well suited for rotating models, since it is com- 
pletely symmetric with respect to J. Therefore we have chosen 
for an alternative limit with the same zeroth and first deriva- 
tive as El at the point of the circular orbits with r — r max 
and with positive J, but with a second derivative which can 
be chosen freely. It is defined by 



E L , (E,J) =E-E, 



c. m a x 



+ 



J, 



c. m a x 



yj <Jc,max^) 



- P(J-Jc, max ) 2 > 0, (A2) 
with the binding energy for circular orbits with positive J at 



Ec max — Vof^max) + n Tmax T (^raai)i 

2 ar 



(A3) 



and the angular momentum at the same point 



Jo. 



dV 

T"max , \ f max ) 

ar 



(A4) 



The parameter /3 is adjustable. For large r, this new limit 
lies closely to circular orbits with positive r, and is much better 
suited for rotating models. I n or der t o av oid unwanted orbits 
with r+ > r max , both limits (|A2| ) and ([A3]) should be combined 
with the extra condition (see also figurefiq) 



E ^ E c ,max • 

The unperturbed models are defined as 



(A5) 



fi = E l , .q x ( 3.84 x 10~ 3 . E A 



+2.22 x 10" 
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(A6) 



/n = -El, o.oo7 x ( 



3.68 x 10" 
+3.15 x 10" 



£ 4 

E 4 pw(J, 2) 



), (A7) 



This article was processed by the author using Springer- Verlag IWgX 
A&A style file 1990. 



